function [f_val_stack, Wald_store, theta_mat]=Wald_CS_calc(alpha,data,theta_grid_array,Farray,giveMoments,giveWeight,giveVarianceEstim)
%Calculate Wald confidence sets. theta_grid_array
%is a cell array, where cell i contains the grid of values we want to
%consider to parameter i.  Here I assume we are interested in confidence
%sets for linear combinations of parameters, and the linear combinations of
%interest are given the Farray.  GiveMoments is used to pass the moment and
%derivate functions, while giveWeight is used to pass the desired weighting
%matrix and giveVarianceEstim passes a function used to calculate
%covariance matrices. alpha sets coverage, and a_vec sets the weight to be used
%in linear combination tests (a_vec should be of the same length as the
%third dimension of Farray).

%Calculate number of parameters
m=length(theta_grid_array);

%Generate vectors consisting of all parameter combinations consistent with
%the parameter grids
theta_mat=theta_grid_array(1);
theta_mat=cell2mat(theta_mat)';

if m>1
    for n=2:m
        theta_old=theta_mat;
        theta_new=theta_grid_array(n);
        theta_new=cell2mat(theta_new)';
        theta_mat=[kron(ones(length(theta_new),1),theta_old) kron(theta_new,ones(length(theta_old),1))];
    end
end

%Define GMM objective
moments=str2func(giveMoments);
W_func=str2func(giveWeight);
g=@(theta) mean(moments(data,theta,0))';
W=@(theta) W_func(data,theta);
obj=@(theta) g(theta)'*W(theta)*g(theta);

%Evaluate GMM objective on parameter space grid to get good starting value
obj_store=zeros(size(theta_mat,1),1);
for n=1:size(theta_mat,1)
    theta_temp=theta_mat(n,:)';
    obj_store(n)=obj(theta_temp);
end
theta_min=theta_mat(obj_store==min(obj_store),:);
theta_min=theta_min(1,:)';

%Calculate GMM estimator
[theta_hat, Sigma_theta_hat] = GMM_Estimator(data,theta_min,giveMoments,giveWeight,giveVarianceEstim);

%Calculate Wald confidence sets for parameters of interest
Wald_store=zeros(size(theta_mat,1),length(Farray));
f_val_stack=cell(size(Farray,3),1);
for k=1:size(Farray,3)
    F=Farray(:,:,k);
    F(:,sum(abs(F),1)==0)=[];
    p=size(F,2);
    Wald=@(theta) size(data,1)*(theta_hat-theta)'*F*(F'*Sigma_theta_hat*F)^-1*F'*(theta_hat-theta);
    for n=1:size(theta_mat,1)
        theta_temp=theta_mat(n,:)';
        Wald_store(n,k)=Wald(theta_temp);
    end
    init_Wald_CS=theta_mat(Wald_store(:,k)<chi2inv(1-alpha,p),:);
    f_CS=init_Wald_CS*F;
    f_val_stack{k}=f_CS;
end

end

